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Abstract 

We study the two-scale asymptotics for a charged beam under the action of a rapidly oscil- 
lating external electric field. After proving the convergence to the correct asymptotic state, we 
develop a numerical method for solving the limit model involving two time scales and validate 
its efficiency for the simulation of long time beam evolution. 
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1 Introduction 

Charged particle beams, as used in beam physics for many applications ranging from heavy ion 
fusion to cancer therapy, often need to be transported on very long distances in particle accelerators. 
Moreover, as particles of the same charge repulse each other, these beams need to be focused using 
external electric or magnetic fields (see, for example, the paper of Filbet and Sonnendriicker [12j for 
a mathematical description of particle beam modelling and numerical simulation). Beam focusing 
can be performed either in linear or circular accelerators using a periodic focusing lattice. 

The aim of this paper is to study the evolution of a beam over a large number of periods. The 
motion of electrically charged particles is governed by several phenomena: in priniis, it is necessary 
to take into account the interactions between the electric fields generated by the particles themselves 
and by the external focusing electromagnetic field. The kinetic approach is a successful method 
suitable to model a system composed by a great number of particles; in our case the modelling 
framework will be given by coupling a kinetic equation with the Maxwell equation. When it is 
possible to disregard the collisions between particles, the kinetic part of the modelling analysis is 
performed by means of the Vlasov collisionless equation [13j . In this paper we will consider only 
non-relativistic long and thin beams; therefore, instead of studying the phenomenon by means 
of the full Vlasov-Maxwell system, we will consider its paraxial approximation. This simplified 
model of the Vlasov-Maxwell equations is particularly adapted to the study of long and thin beams 
and describes the evolution of a charged particle beam in the plane transverse to its direction of 
propagation. 

The key assumptions of the paraxial approximation are the following: 
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1. we suppose that there is a predominant length scale, namely the longitudinal length of the 
beam, whereas the transverse thickness of the beam is negligible with respect to its longitu- 
dinal length; 

2. we will assume that the beam has already reached its stationary state. 

Such assumptions are often verified in applications: in a standard accelerator, the orders of mag- 
nitude of the different lengths taken into account verify the first assumption; moreover it is often 
interesting to study stable (in time) configurations of the beam, and therefore the second assump- 
tion applies. In the paraxial framework, the effects of the self consistent magnetic and electric 
fields can be both taken into account by solving a single Poisson equation. Hence the model be- 
comes similar to the Vlasov-Poisson system in the transverse plane with respect to the accelerator 
axis. In the paraxial model, the z coordinate plays the role of time. We will switch here to more 
conventional notations and consider the time-dependent two-dimensional Vlasov-Poisson system. 
We assume moreover, in accordance with most experiments, that the average external field is large 
compared to the self-consistent field. More details about the derivation of this model can be found 
in |1] or [12], whereas a more complete description on accelerators' physics can be found in the 
book by Davidson and Qin [3]. 

In the framework we just precised, the Vlasov-Poisson system takes the form 

(9f 1 

^ + -V • V,./, + (E, + H^) • V„/, = 0, 
at e 



= -p£(f,x), Pe{t,X.)= /e(t,X,v)dv, 



_ /,(t = 0,x,v) = /o, 

where fg = /e(t,x, v) is the distribution function, t £ [0,T) for some T < cxd, x = (xi,X2) € is 
the position vector and v = {vi,V2) G is the velocity vector. We will denote by O = x M^. 

As said before, we are working in a stationary setting. In this context, the variable t does not 
represent, from a physical point of view, a time variable, but rather the longitudinal coordinate. 
We have nevertheless chosen to use the variable t because of the great similarity of the paraxial 
approximation with respect to a two-dimensional (in space) time-dependent Vlasov equation. The 
electric field H'^(t,x) is supposed to be external, and has the form 

3^(t,x) = - ^Hox + Hi(^uji^^ Kx, 

where Hq is a positive constant and Hi is 27r— periodic, with mean value 0. K is a 2 x 2 symmetrical 
isometric matrix. In the two most common physical situations, K = /orK=(^ The 



.0 -1, 

parameter e is a scaling parameter acting on the "time" scale (i.e. the longitudinal scale). Its 
physical meaning will be detailed in the next section. 

The form of the external electric field will satisfy some relevant physical requirements. The field 
which models the periodic focusing lattice needs indeed to be linear in x to compensate for the self 
field which is also linear in x for a uniform beam but with the opposite sign. The variation of the 
field with respect to t is assumed to be a small amplitude field oscillating at the same period as 
the lattice. 
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Under additional physical hypotheses, a simplified version of ([T]) may be considered. These 
assumptions consist, in the case when K = /, in considering that the beam is axisymmetric, i.e. 
invariant under rotation in M^. Then, by writing ([Ha) in polar coordinates (r, 6), such that r = |x|, 
xi = rcos{6) and X2 = rsin(^) we get 



dfe i dfe , ^ 



dt 



dr 



dVr 



0, 



where Vr is the projection of the velocity on radius x. In other words, Vr = v-x/r and vq = v-x-*"/?" 
where x-*- = {—X2,xi). If finally, is assumed to be concentrated in angular momentum rvo = 
(this property is achieved if the initial condition /o is concentrated in angular momentum), system 
^ yields 

19/.,^ , „ dfe 



0, 



1 d{r^r 



Peit,r), pe{t,r)= / fe{t,r,Vr)du 



(3) 



r dr 

fe{t = ^,r,Vr) = fo, 

where fe = fe{t,r,Vr) for t G [0, T), r G and Vr G M. The external force writes here 



S/(t,r) = {-^Ho + Hi(^LJi^-^ 



System ([3]) is naturally set for r G M+. Nevertheless, we can consider that r G M using the 
convention fe{t,r,Vr) = fe{t, —r, —Vr) and Sr^{t,r) = Hr'^(t, — r). This is what we do in the 
following. 

The main purpose of the paper is the study of systems ^ and ([3]) in the limit as e ^ 0. We 
have chosen to work in the framework of the so-called two-scale convergence, an homogenization 
technique introduced by N'guetseng [M] and subsequently developed by Allaire [2] and Frenod, 
Raviart and Sonnendriicker [9], which essentially says that /e(t,x, v) is close to F{t,t/e,x,v) for 
a given profile F. This notion is stronger than the usual weak or weak* convergence and is very 
adapted to the study of asymptotic behaviour of transport equations with rapidly oscillating terms, 
as was shown in O [HI [9l [lOl [11]. On the other hand, we develop a numerical method in order to 
compute an approximation of fe- This method consists, in the spirit suggested in Frenod, Raviart 
and Sonnendriicker [9j and Frenod [6] and explored in Ailliot, Frenod and Monbet [1] and in Frenod, 
Mouton and Sonnendriicker [3, in computing the profile F by discretizing the equation it satisfies. 
Then an approximation of fe is reconstructed using fe{t,x,v) F{t,t/e,x,v). For simplicity 
reasons, the considered method is implemented in the simplified case of system ([3|). Nevertheless, 
it is easy to see that it may be adapted in a straightforward way to system ([1]) and others situations 
where multiple time scales are involved. 

The paper is organized as follows: in Section [2] we obtain the scaled Vlasov-Poisson system ([T|) 
from the standard Vlasov-Poisson system by means of a scaling procedure. Section[3|is then devoted 
to the derivation of some useful properties of the solution of ([T|), whereas Section [5| is dedicated to 
homogenization of systems ([T|) and ([3|) by means of the two-scale convergence. Finally, in Section [5] 
we will implement and test the previously described numerical method based on the homogenization 
analysis. 
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2 Scaling of the Vlasov equation: the paraxial approximation 



We present in this paragraph the scahng procedure leading to Equation ([T]). We first note that, 
in the formulation of the problem, there are different scales: the length of the accelerator L, the 
transverse radius of the beam A and the period of the oscillating external field /. In a standard 
accelerator, the orders of magnitude of such quantities are: L ~ 10^ m; / ~ 10~^ m and, finally, 
A ~ 10~^ m. The standard three-dimensional Vlasov-Poisson system, in presence of an external 
electric field H, is given by 



^ + V • V,/ + ^ (E + H) • V,/ 

ot m 



0, 



E = 



-V$, 



1 



— p(i,x) 
£0 



p(x,i) = q 



(4) 



/(t = 0,x,v) =/o, 



where / = /(t,x, v) is the distribution function, t e [0, T) for some T < cxd, x = {xi,X2,X3) S 
is the position vector and v = (f i, ^2, "^3) G is the velocity vector. 

The external electric field H = H(x) is supposed to be independent of time and periodic with 
respect to the variable X3. Under the paraxial approximation, / is supposed to be stationary 
with respect to time and the velocity of the particles with respect to the longitudinal direction 
of the beam (which will be 63 throughout the whole paper) is constant, that is = Vb- These 
assumptions will enable us to eliminate the time derivative in the first equation of (jlj) and suppose 
that /(t,x, v) = f{x,vi,V2) S{vb — V3). In order to obtain the dimensionless version of the 
Vlasov-Poisson system we define the new variables x'^, v[ by: 



Ax-, 



Vb Vi, 



1, 2 



X3 = LX3, V3 = Vb. 

Moreover, the dimensionless density /', external electric field 3', self-consistent electric field E' and 
its potential are defined by 



E = EE\ 



<|) = (J) <|)'. 



Under such a change of variables, we obtain that the Vlasov-Poisson system @ can be written in 
the following form: 



fvbdf^ fvb^ ,df' qEf^f df 



L dx'n 



EE' 



A 

A a^' 



mvb 



4=1 



0, 



1, 2: 



EE'. 



L dx'^ ' 



(5) 



A2 ^ V 9af 

i=l ' 



I2 &f 



+ 



-P'(x'), 



P ix 



£0 



/ /'(x',v')dv' 



ff'{t = 0,x',v') = fo{Xx'i,Xx'2,Lx'3,VbVl,VbV2). 
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It is natural to choose, as unit measure of the electric field 



E 



mv 



b . 



qL 



moreover, since the third component of the electric field will play no role, its unity measure induces 
that 



<^ = EX 



qL 



Finally, if we assume to have a small self consistent electric field, we are forced to work with a 
small density of particles. Therefore, we will assume that the dimensionless form of the density / 
is given by / = / /', with 

/ 



q^XL' 



As mentioned in the introduction, we consider here low intensity beams so that the external 
electric field is much stronger than the self-consistent field. We translate this, introducing a small 
parameter e, by saying that r./E = 1/e and assuming that H' writes 



where Hi is 27r— periodic and for a 2 x 2 isometric matrix K 



Gi 




and G2 




On the other hand, we also suppose that the ratio between the average transverse radius of the 
beam and the longitudinal length of the accelerator A/L = e, and, that the ratio between the period 
of the external electric field and the macroscopic quantity l/L = e. Forgetting the terms in e^, we 
have then obtained the dimensionless Vlasov-Poisson system in the paraxial approximation: 



E[ = -d^'/dx[, 



. , dx: 



1,2, 



P'(x') 



(6) 



/'(x',v')dv'. 



/'(t = 0,x',v') = /,^, 



1 



with 



Denoting the variable X3 with t (in order to have in ^ only two-dimensional vectors) and elimi- 
nating all the primed variables, we finally obtain the scaled Vlasov-Poisson paraxial system ([1]). 
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3 Main properties of the solution 

We recall that, from now on, all the vectors are intended to be two-dimensional. 

In the sequel, we will use heavily the property of boundedness of the solution with respect to 
some norms. This is guaranteed by the following lemma: 

Lemma 3.1 Let fi, he the solution of the Vlasov-Poisson system ([I]), with non-negative (a.e.) 
initial data fo of class {L^ H L^){^1), p > 2 and such that the moment of order two is finite: 

/ (|x|^ + |v|^)/odx(iv < +00. 
Jq 

Then f^ is bounded in L°°{[0,T]; L'P{il.)) uniformly with respect to t. Moreover 

||(|x|2/e + |v|V.)llL-([0,T];Li(n)) <Ce"^, (7) 



with 



Finally 



^ll-f^llloo 

a 



min{l , Hq}' 



\pe{x,t)\\ ^ , 3_,, < Ce°^/3 (8) 



for some constant C. 

Proof: We multiply the Vlasov equation ([1^) by fi~^ and integrate in x and v. We obtain 

ll/'^llL°°([0,T];LP(n)) < C, 

for some constant C. This proves the first part of the lemma. 

In order to prove the other statements, we multiply the Vlasov equation ([1^) by |vp, and 
integrate with respect to x and v. We get 



dt 



[ /e|v|2 dxdv-2 [ J • (Ee + H^) dx = 0, 



(9) 



where 



J(x,i) = / v/,dv. 
Now, integrating the Vlasov equation in v gives the continuity equation 

^ + -V-J = 0. (10) 

at £ 

Thanks to the continuity equation pU|). we note that 

3-Eedx = - [ 3-V<^,dx= [ V-J$edx = -e/ ^^edx. (11) 
Using now the Poisson equation, we get 

-4/ {V^efdx = -[ 4- (A$,) dx = / ^^edx. (12) 
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Coupling pT]) and ([T2|) we deduce: 

-2 / J-Ee(ix = e^ / (V$e)2dx. 

Jm2 at J^2 

On the other hand, the external electric field H"^ can be also derived from the potential 

P(x) = ^i7o|x|2-i//i(^cJi^^ Kx-x. 

A variant of the procedure used for the self-consistent field allows us to deduce 

dp, 



(13) 



J • dx = 2 / J • VP dx 



V • J P dx = 2e 



dt 



Ho [ |xp^dx-e / i^if wi- ) Kx -x^dx, 

Jr? dt \ ej dt 



and therefore (191) becomes 



/ fe\v\^dvdx + e / (V$e)^ dx + / |x|VeC^x 
J J JRl 



eHi uji 



t 



Kx • X— ^ dx. 

dt 

If we multiply the Vlasov Equation ([1^) by |xp and then integrate with respect to x and v, we 
deduce that 



Kx • X — ^ dx 

dt 



/ Kx • X V • Vxfe dx dv = 2 / (Kx • v) dx dv 



and thanks to the elementary inequality 2Kx • v < (|Kxp + |vp) < (|xp + |vp), we obtain that 



d_ 

dt 



/ /elvl^dvdx + e / (V^>e)^dx + Fo / l^l^Pedx 



< 



< e\\H- 



1 ||oo 



|xp/£dxdv+ / Ivp/gdxdv 



Finally 
with 



|x|V. + |v|V.)llL-([0,T];Li(n)) <Ce"^ 
^||-f^l||oo 



a 



min{l , Ho} 

which proves the second statement of the lemma. 

The proof of the bound on /?£ is a straightforward consequence of the following classical estimate 
(see for instance [IT] (lemma 4.4)): 



[ \peix,t)f/^dx<cJ [ Uefdxdv)^ ([ \v\^fedxdv 



1/2 



Hence the proof of the lemma is complete. 



□ 
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4 Homogenization of the paraxial approximation 



Our goal consists in deducing the equations satisfied by the hmit of as e — > 0. The homogeniza- 
tion of partial differential equations with rapidly periodic coefficients, as in the present case, can 
be fruitfully studied applying the two-scale convergence, introduced by N'guetseng [Hj and subse- 
quently developed by Allaire [2j and Frenod, Raviart and Sonnendriicker [9j and which essentially 
stands in the following theorem. 

Theorem 4.1 Let Q be a regular subset o/M", X a Banach space and X' its dual space, with { , ) 
as duality bracket. Let 1 < p < +oo and p' be such that ^ + ^ = 1 and let Y = [0,ai] x ■ ■ ■ x [0, a„] 
defined for finite real numbers ai . 

If a sequence (n^) = (us{q)) is bounded in L'^ {Q;X') then there exists a function uq = UQ{q,y) G 
LP {QxY;X') such that, up to a subsequence, 

lim^ I) )dq = j ^ j^{uo{q,y),v {q,y)) dqdy, 

for any function u £ LP{Q; Cy(M";X)), where Cy(M";X) stands for the space of continuous 
Y— periodic functions on M" with values in X. We say that (^^(x)) two-scale converges to uo{x,y). 

More precisely, we will use a theorem of Frenod and Sonnendriicker [TT] , which gives a generic 
framework, in the context of two-scale convergence, for linearly perturbed conservation laws of the 
form 

— 3^ + • V^Ue + -L • V^Ue = 0, 
ot e 



(14) 

. Ueit = 0) = Uq. 

In this system, u'^ = u'^{t,x), t £ [0,T) for some T < oo and x G M" = O. Moreover it is assumed 
that, for all e > 0, Vc • A*" = 0, and that, for some q > 1, A^ = A^(t,x) two-scale converges to 
A = A{t, T, x) E L~([0, T] X [0, 0] ; W^''}{K)) for all compact sets K G M". Finally, L = Mx where 
M is a real nxn matrix with constant entries, satisfying trM = and such that e"^^ is 0— periodic. 
Under these conditions, the following theorem is proved in |llj : 

Theorem 4.2 Under the assumptions above, if moreover the sequence (u^) of solution of ( | j^P 
satisfies 

\\Ue\\L°°i[(},T];LP(0)) < C, (15) 

for some p > 1 such that 1/p + 1/q' < 1, where 1/q' = max{l/g — 1/n, 0}, then, extracting a 
subsequence, Us two-scale converges to a profile 

U G L°°([0,T] X [0,9];LP{O)). 

Moreover we have 

U{t,T,x) = Uo{t,e-^^'x), (16) 

where Uq = UQ{t^y) is solution to 

dUo 



dt 



^ + / e-'^^Ait, a, e"*^y) da ■ VyUo = 0, 

^ (17) 



Uoit = 0) = ^no. 
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The problem of homogenizing system ([T]) enters the framework we just presented with x replaced 
with (x,v) E ri, with 

,E,(t,x)+Fi(cui^)Kx^ 



(18) 



and with 



It is an easy game to see that 

where 
and 

We can now prove: 



or in other words M 



/ 

-Hoi 



7^l(T) 7^2(T)/u;o 
-wo7^2(r) 7^l(T) 



with LJQ = y' Hq, 



(19) 



n2{s) 



COS {ijJqs) 

COS {ijJqS) 

sin {ijJqs) 

sin {ijJos) 



Theorem 4.3 Let /q satisfy the hypotheses of Theorem Vj.l\ Then, if we consider a sequence of 
solutions (/e,Ee) depending on e, extracting a subsequence, we have that f^ two-scale converges to 
FGL°°([0,r] X [0,^];L2(0)) and^e two-scale converges to^GL°°([0,r] x [Q,'^]-W^^^'^{Rl)). 



Moreover, there exists a function G = G(i,y,u) E L°°([0,T];L (0)) such that 

F{t,T,x,v) = G (^^,7^l(-r)x + -^7^2(-T)v, -(Jo^2(-t)x + 7^l(-T)v 

The pair (G, 8) is finally the solution of 

' dG 1 r'^'^/'^o 
-K7 + — n2{-a)S{t,a,niia)y + —n2ia)u)da-VyG 

Ot LOo Jo Wo 

r-2-K/ujo 

+ / 7^l(-fT)f(^,a,7^l((T)y + — 7^2(cT)u)(i^7•V„G = 0, 
Jo '^0 

G{t = 0) = ^/o, 

if ui/ujo ^ Q or if Hi = and 
' dG 1 /■2'^/'^o 



(20) 



1 



+ — / ^2(-fT) f(^,a,7^l((7)y + — 7^2(^7)u) 
dt cjo Jo L ^0 

+ ^/^l(L^lC7)K(7^l((T)y + — 7^2(a)u) 
zvr Wo ■ 

+ / 7^l(-^7) ^(^,a,7^l(a)y + — 7^2(^T)u) 

+ ^//l(L^lC7)K(7^l(a)y + — 7^2('T)u) 
27r Wo ■ 

G(t = 0) = ^/o, 



da-VyG 



da-VuG = 0, 



(21) 
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if LOi/ujQ G N, where E = £{t,T,x) is solution to 
( S = -V^, 

r 1 (22) 

-A^- = / 7^l(-r)x + — 7^2(-T)v , -a;o^2(-r)x + 7^l(-r)v) dv. 

Proof: The deduction of this theorem follows directly from theorem 14.21 using the expression (jl9p 
of e'^''^ once the two-scale convergence in A'' given by (llSh and in the Poisson equation ([it) is 
achieved. 

As a direct consequence of estimates ([7]) and ([8]) and of the regularization properties of the 
Laplace operator, we deduce that is bounded in L°°([0, T]; Ty^'^/^(M^)), then, extracting a 
subsequence, two-scale converges to £: = -V^- G L°°([0,r] x [0, W^1'3/2(r2)). 
In order to pass to the two-scale limit in the Poisson equation, we multiply it by a test function 
v{t, |,x) such that r ^ ^{t^T,^.) is 27r/LiJo~periodic, to give 

nV<i>e{t,x) ■Vu{t,-,x)dxdt = I I fe{t,x,v)u{t,-,x)dxdvdt, 
_ _;2 e Jo Jn e 

in which we pass to the limit to obtain 

/"°/ / V'^{t,x) ■Viy{t,T,x)dxdtdT = /"°/ / F{t,T,x,w)u{t,T,x) dxd^dtdr 
Jo Jo Jr^ Jo Jo Jn 

= I / G (t, 7?.i(— r)x + 7^2(— ''")v/ciJo , — i^o^2(— ''")x + 7^i(— r)v) z^(t, r, x) dxdvcitdT, 
Jo Jo Jq 

which is the weak formulation of (j22p . 

On the other hand, if ui/ujo ^ Q, for any regular function which is 27r/ci;o— periodic in r, the 
product 

-)Kx • i^{t, -, 

since r i-^ Hi{ujit) is 27r/a;i— periodic with mean value 0. This yields equation ([20]) . 
If ioi/ujo G N, then r i— > Hi{loit) is 27r/a;o— periodic and 

n//i(a;i-)Kx • -, x) (ixdt ^ / °/ / Hi{ujit)}^x ■ u{t,T,x) dxdtdr, 
yielding (j2ip and then ending the proof. □ 



Hi{uJi-)'K.x ■ u{t, -, x) dxfit —>■ 0, 



Remark 4.4 I^e can bring back the case loi/ujq S Q \ N to the case loi/ujq S N by finding two 
integers k and I such that uj'q = uJo/k and uj[ = uJi/l satisfies uj'i/uj'q £ N and replacing in ()^ip 

/ da by da. 
Jo Jo 

We now turn to homogenizing system ([3]). For simplicity, we assume here that loq = Hq = 1. 
We also assume that /o satisfies 

/o > 0,/o G (L^ n LP) (M2;rdrd^;^) for j5 > 2 and / {r^ + v'^)fordrdvr < +oo. (23) 

JR2 
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The considered system enters the framework presented in the beginning of the section with x 
replaced by {r,Vr) G with 



M 





Ere{t,r)+Hi{LUil) rj ' 

tM 



Vr 

— r 



1\ M ^ ( cos(r) sin(r) 

-1 Oy ' V -sin(r) cos(t) J ' 

and the solution satisfies the needed estimates to deduce the following Theorem. 



(24) 
(25) 



Theorem 4.5 Under the assumptions above, extracting from a sequence of solution (/gjE^.^) to 
(0) a subsequence, we deduce that fi. two-scale converges to F £ L°°([0,T] x [0,2tt]; L'^ (M?;rdrdvr)) 
and Ere two-scale converges to £ £ L°°{[0,T] x [0, 27r]; ^/^'^/^(M; rdr)). 
Moreover, there exists a function G = G{t,q,Ur) G L'^ {[0,T]; L'^ [M.'^ ; qdqdur)) such that 



F{t, T, r, Vr) = G{t, cos(T)r — sin(r)?;r, sin(T)r + cos(r)?;r), 

and G is solution to 

( dG 



(26) 



dt 



»27r QQ 

+ / —sm.{a)£r{t,a,cos{a)q + svii{a)ur)da 







+ 



oq 

I cos{a) £r{t-,(y-,cos{a)q + sm{a)ur) da — 
Jo , 



dG 

dUr 



0, 



(27) 



ifuJi/uJo 



G{t = 0) = -/o, 



or if Hi = and 

2tt 



dG f ( 

— — h / — sin(cT) I (?r(i, o", cos(cr)g + sin((T)nr) 
ot Jo V 

+ —Hi{uJia)(^cos{a)q-\-sm{a)ur)^ da 

+ J cos{a)^£r{t,a,cos{a)q -\- sm{a)ur) 

1 \ dG 

+ —Hi{u!ia){cos{a)q + sin(cr)ur) ) da - — = 0, 

ZTT J OUr 



dG 
dq 



(28) 



^ = 0) = ^/O' 
ifuJi/ujQ G N where £r = £r{t,T,r,Vr) is given by 
1 d{r£r) 



dr 



T{t,T,r) = / G[t,cos{T)r — sm{T)vr,sm{T)r -\- cos{T)vr) dvr 
Jr 



(29) 



5 The two scale PIC solver 



In this section, we develop a two-scale PIC-method, tailored to approximate fs very efficiently on 
long time scales, in the simplified case of axisymmetric beams. The strategy consists in computing 
the solution F to (l26])-(l27])-(l29]) or ([26])- ([28])- ([29]) and then to approach the solution fe of © by 
F(t, i/e, x, v). The advantage of proceeding in such a way is that the solution of (f26]) - (f27|) - (f29]) or 
(j26p -()28p- ()29p does not contain 1/e— frequency oscillations. As a consequence, a much larger time 
step may be used in the numerical method. First, we present the implemented algorithm. Then, 
we compare the solution obtained with our method with the solution directly computed from 
([3]). Finally, we compare the performances of both methods. 
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5.1 Description of the numerical method 

The Particle In Cell (PIC) method for the Vlasov-Poisson (or Vlasov-Maxwell) equations consists 
in approximating the distribution function defined in phase space by a meshless particle method 
and the electric field on a grid of the physical space only. In addition, adequate interpolation and 
charge deposition methods are used to transfer the needed quantities between grid and particles. 

For the two-scale method we want to develop, the Vlasov-Poisson equation are replaced by 
equations ([26])- ([271)- ([29]) or (f26D -(l28D-(l29D. 

The key point of the algorithm is the computation of function G solution to ([27|) - ([29]) or ([28|) - 
(j29l) at time = ti + At knowing it at time ti. We give the details of the method only in the 
case of system (p7|) - ([29|) . knowing that in the case of system ([28|) - (p9]) ad-hoc terms need to be 
managed. 

As usual in a PIC-method, G is approximated by the following Dirac mass sum 

N 

GNiq, u,t) = Y^ Wk5{q - Qk{t))S{u - Uk{t)), (30) 
k=l 

where (Qfc(t), Uk{t)) is the position in phase space of macro-particle k which moves along a charac- 
teristic curve of the first order PDE (|27p . Hence the job is reduced to compute the macro-particle 
positions {Q^j^^ , ^k^^) time = ti + /S.t from their positions (Q^, U]^) at time tu knowing they 
are solutions to 



dQk _ _ 
dt ~ 

dm '■2- 



dt 



/ sin{a)£r{t,a,cos{a)Qk + sm{a)Uk)da, Qkik) = Q\, (31) 

Jo 

1-2-K 

/ cos((t) <£",.(*, cj, cos(cr)(5fc + sin(cr)[/fc)d(T, Uk{ti) = Ul- (32) 

Jo 



These characteristic equations are a lot more complex to deal with than those of the standard 
Vlasov-Poisson system which read ^ ~ ^ -> ^ ~ E{X,t). Indeed, in our case Q and U are 
coupled in each equation and an integral term needs to be approximated. 

Because of the form of the right hand side in (j3ip and (j32p all along the algorithm, we need to 
compute values of the two-scale electric field £r generated by a given macro-particle distribution 
{Qk,Uk)k=i,..,N- The tedious step while computing £r in grid point {qi)i=i,..,A is the computation 
of the right hand side T{t, am, Qi) of (p9]) . Indeed, T(t,am,qi) involves the integral of the particle 
distribution on the oblique line which is the range of the vertical line [q = qi] by the rotation 
g-(T„A/ (igfined by (j25p . Hence we have to apply this rotation to each line [q = qi] and to project 
the particles (Q^, Uj^)k=i,..,N on the resulting oblique lines. Summing then the projection result on 
the oblique line associated with qi yields the value of T(t;, dm, qi)- Another way to obtain this value 
consists in applying the rotation e'^™-^^ to the particles, then projecting this rotation result on lines 
[q = qi] and summing. Once T(t;, cjm, qi) is known in each qi, the computation of the £r{ti, am, qi) 
are straightforward using any classical Poisson numerical solver. 

The first step of the computation of {Q^^^ ,U^]^^) consists in replacing the integrals above by 
p— node quadrature formula. As we approximate the integral of a periodic function over one period, 
the trapezoidal rule is optimal and will yield very accurate results for as few quadrature points as 
are needed to resolve the oscillations of the function. 

Then, the equations for {Qk, Uk) become 

dQk ^ 

-j^ = -'^^mSm{am)£rit,am,cos{am)Qk + sm{am)Uk), Qk{ti) = Q[, (33) 

m=l 
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dUk 
dt 



jmCos{am) £r{t, crm,cos{am)Qk + sm{am)Uk), UkiU) = Ul- 



(34) 



m=l 



Then, we solve (j33|l - (pl|l using the classical Runge-Kutta method: 










1/2 


1/2 




1/2 





1/2 


1 





1 




1/6 


1/3 1/3 1/6 



(35) 



which gives the following scheme when applied to the computation of the approximation y'+^ of 
the value of y solution to dy/dt = K{t, y) at time ti + At knowing its approximation at time ti: 



ti,i=ti, y^^^=y^ 

ti,2 = *i + y''^ = y^ + \l\ with /I = At K(t/,i, y 



At 

ti,3 = ti + —, y- 



1,3 



1 



ti^i = ti + At, y''^ = y' + /^ with = At K{ti^3, y''^), 



' + with/2 = Ati^(t;,2,y''') 



(36) 



y'^' = y' + \i' + \i' + \i' + with = Atm^,J'% 

6 3 3 d 

Applying this scheme to our problem consists in replacing in the formula above y by (Qfc, Uk) and 
computing K using the result of a Poisson solver. In other words, we have to compute as 
follows: 

1^2 _ /^i 1 



gr = Qi + -/\ with 



= Ati - ^ 7mSin((Tm,)<?r(i/,0-m,COs((Tm)Qi + s\n{a.m)u[] 



tL2 



(37) 



and something similar for C/^ . 

In order to achieve this, we compute the value of £r in (t;, o"m, cos((Tm)<5fc + sin((Tm)f/^) by 
interpolating the value of £riti,crm, Qi) known on the grid {qi)i=i..,A as soon as it has been computed 
solving the Poisson equation (p9|) associated with the particle distribution (Q^,, Ul) by the procedure 
described above. 

The following step of the Runge-Kutta method consists in computing Qjl defined by 



= Qi + with 



At 



7,2x 

k J/' 



(38) 



m=l 
?2^^.. I Af 



where the value of £'^{ti + ^ , Cm, cos{am)Q^jf + sin{am)U^k'^) is obtained by interpolation of £'^(t/ + 

At 
2 



r«,2N 



'-,am,qi) which is computed as previously from the {Q^'jf ,Uj^'^)k=i,..,N particle distribution. 
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Then we compute 



Qi;^ = Qi + /^ with 




p 



(39) 



where £^{t + ■^) is computed from particle positions {Q^if, U\f')k=i,..,N ■ 
Finally, Q^^^ is obtained by the following formula: 




(40) 



m=l 



where I^, and are defined above and where £ji{ti + At) is computed from particle positions 



5.2 Validation of the two-scale solver 

The linear case. First in order to check our implementation we consider the case when the self- 
consistent electric field vanishes. Then by choosing adequately the function Hi we can compute an 
analytical solution of the two-scale model that we can compare with our code and with the solution 
given by the usual PIC solver. 

Let us first consider the non resonant case by choosing Hi{t) = cos(4-v/2t). In this case, as 
the self-consistent electric field vanishes, G is stationary and the beam will only move through the 
rotation transforming G to F. We check on figure[T]that this is indeed the case. Let us now consider 
the resonant case by choosing Hi{t) = cos^(nt) with n > 2. Then a straightforward computation 
yields the equation satisfied by G 



which yields an explicit solution of G and also an explicit solution of F using formula (j26p . Com- 
paring this exact solution with the solution computed by our solver we can check its accuracy. Let 
us first mention, that our quadrature formula yields the exact result up to machine accuracy for any 
number of quadrature points greater or equal to seven for n = 2 in the definition of Hi . Of course 
when the chosen n is larger the function oscillates more and more quadrature points are needed. 
For n = 7, for example, the exact result, up to machine accuracy, is obtained with 17 quadrature 
points. We can thus conclude that as long as the oscillations are resolved our quadrature rule is 
very accurate. We also check that our RK solver is of order 4 in At as expected. In figure [21 we 
check the global behavior for a whole beam in comparison to the usual PIC solver and the results 
are very satisfying. The phase shift appearing in the two-scale limit is indeed also observed in the 
usual PIC code. Note that because the usual PIC code needs to resolve the fast time scale, the 
difference in time step between of the two solvers is of the order of e, so that the efficiency of the 
two-scale PIC solver compared to the usual PIC solver is better for small values of e. 




The characteristics of this equation can be computed explicitly: 



Q{t) = Qo cos - - C/q sin -, U (t) = Qq sin - + Uo cos - 
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Figure 1: Beam simulation without self-consistent electric field in the non resonant case {Hi{t) = 
cos(4^/2i)) with a usual PIC method (left) and a two-scale PIC method (right) for e = 0.01 at time 
6.28. 



.3J.H I |.| mi In ij 1 tij I fi III It 



Phase space at time = 6.2B 
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Figure 2: Beam simulation without self-consistent electric field in the resonant case (Hi{t) = 
cos^(2t)) with a usual PIC method (left) and a two-scale PIC method (right) for e = 0.01 at time 
6.28. 
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Phase space at time = . 031 Phase space at time =34 . 558 Phase space at time = 69 . 115 
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Figure 3: Beam simulation with a usual PIC method and a two-scale PIC method for e = 0.01. Left: 
beam at time 0.031, center: beam at time 34.558, right : beam at time 69.115. Top : Simulation 
provided with the usual PIC method, bottom: Simulation provided with the two-scale PIC method. 



The non linear case. In this case, as no analytical solution is available, we shall compare our 
solver with a traditional PIC solver resolving the small time scale. 

Let us first consider a case of a fairly small e on which we test our two-scale PIC method, where 
the two scale solution should be fairly close to the solution given by the traditional PIC solver. We 
take £ = 0.01, and choose as an initial distribution 

h{r,Vr) = exp( - ^^X[-f).7bfl.7S\{'r-), (41) 

with thermal velocity Vth = 0.0727518214392 and where X[-o.75,o.75]('^) = 1 if r G [-0.75,0.75] 
and otherwise and in considering function Hi as being identically zero. This corresponds to a 
semi-Gaussian beam often used in accelerator physics. We consider in all our simulations beams 
which are willingly unmatched in order to assess to what extent our two-scale PIC solver can follow 
the complex structures that are developing. 

We use a 15-node composed trapezoidal quadrature formula. The results are given in figure [3] 
where the horizontal axis is the r axis and the vertical axis the Vr axis. The top line shows the time 
evolution of the beam simulated with a usual standard PIC method and the bottom line shows the 
same simulation with the two-scale PIC method just built. We can see that the two simulations 
coincide with a high degree of accuracy. The time step in the two-scale PIC method is e times 
larger as in the usual PIC method. The simulations were both made with an Intel Core 2 Duo 
processor (2.33 GHz) under Mac OS X 10.4.10 (8R2218) system. The needed CPU time for the 
simulation with the usual standard PIC method is 4320.382 seconds and is 549.197 seconds with 
the two-scale PIC method. 
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Figure 4: Beam simulation with an external force with no mean effect at time 9.30. Left: with the 
usual PIC method. Right: with the two-scale PIC method. 
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Figure 5: Beam simulation with an external force with no mean effect at time 9.52. Left: with the 
usual PIC method. Right: with the two-scale PIC method. 

The three next tests are done with e = 0.1 and with initial distribution 

fo{r,Vr) = -j=—exp(^ - ^^^X[~1.83271471003,1.83271471003] (0> (42) 

with thermal velocity Vth = 0.0727518214392. In the case when 

Hi{ujit) = cos{t), (43) 



the mean effect of it is zero in the sense that terms containing Hi appearing in (j28p are both zero. 

Numerical results concerning this case are given in figures [4] and [5l On the left of both figures is 
shown the simulation result with a usual PIC method. On the right is shown the simulation result 
with the two-scale PIC method. We can notice that if the beam configuration are very similar at 
time 9.30 (figure H]) the two-scale PIC method is a bit in advance at time 9.52 (figure E]). This 
illustrates a finite e effect making that the two-scale PIC method is sometimes late and sometimes 
in advance with respect to the usual PIC method. 

The case when 

Hi{uJiT) =cos^{t), (44) 
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Figure 6: Left: beam at time 6.26, center left: beam at time 18.85, center right : beam at time 
31.42, right : beam at time 72.26. Top : Simulation provided with the usual PIC method, bottom: 
Simulation provided with the two-scale PIC method. 



generates a real effect on the long term. In other words, the terms containing Hi in ()28p are 



2tt 



2w 



sin((7) cos^(cj) ( cos{a)q + sm{a)ur) da 



-Ur 



1 f'^'^ 3 
— / cos((7) cos^(cr) ( cos((t)(7 + sin((T)nr.) do" = -g. 



(45) 
(46) 



The simulation results are given in figure El The top line shows the time evolution of the beam 
simulated with a usual standard PIC method and the bottom line shows the same simulation with 
the two-scale PIC method just built. The simulations coincide with a good degree of accuracy also 
in this case. 

Finally we consider the case when 



Hi{ujit) = cos(2r), 
which has a defocusing effect on the beam as 



1 f Ur 

— / sin(cr) cos(2(t) ( cos((t)(7 + sin((T)ur) dfj = 
2vr Jo 4 

1 f'^^ q 

— / cos((t) cos(2(t) ( cos((T)g + sm{a)ur) da = -. 
2vr Jo ^4 



(47) 

(48) 
(49) 



The result shown in figure [7] shows that also in this case the two-scale PIC method gives results 
that are qualitatively and quantitatively very close to those of the usual PIC method. 
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Figure 7: Beam simulation with an external de-focusing force at time 2.20. Left: with the usual 
PIC method. Right: with the two-scale PIC method. 

6 Conclusions and perspectives 

We have developed and validated a new PIC solver that can deal very efficiently with a problem 
involving two time scales. The two-scale solver is based on homogenized equations that have 
been obtained analytically in our application. For small values of e the new solver can follow as 
accurately as the original solver the complex structure of the particle distribution that is generated 
in a mismatched particle beam in an accelerator at a small fraction of the cost of the usual PIC 
solver. These numerical results obtained in ID are very promising and the extension to more 
dimensions as well as the used of semi-Lagrangian or other type of solvers for the Vlasov-like 
equation instead of a PIC solver deserve to be investigated for accelerator physics as well as for 
other applications, e.g. strongly magnetized plasmas where the fast motion along the magnetic field 
lines provides a fast time scale. On the other hand, a similar two time scale method can probably 
also be applied in cases where the homogenized equation can only be computed numerically. 
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